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Abstract 

Transparent boundary conditions for the time-dependent Schrodinger equa¬ 
tion are implemented using the R-matrix method. The employed scattering 
formalism is suitable for describing open quantum systems and provides the 
framework for the time-dependent coherent transport. Transmission and re¬ 
flection of wave functions at the edges of a flnite quantum system are essential 
for an accurate and eflicient description of the time-dependent processes on 
large time scales. We detail the computational method and point out the nu¬ 
merical advantages stemming from the open system approach based on the 
R-matrix formalism. The approach is used here to describe time-dependent 
transport across nanostructured interfaces relevant for photovoltaic applica¬ 
tions. 
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1. Introduction 


By continuing miniaturization the current electronic devices have already 
reached length scales of only a few tens of nanometers and, in the past few 
years, quantum mechanical approaches have been extensively used for the 


modelling of the electron transport down to molecular scale [Ij, l2l. In partic¬ 
ular, there is a lot of interest for an efficient description of time-dependent 
coherent transport and examples of physical systems may include nanoscopic 
antennas j^, electron and hole transport through nanostructured interfaces 
with applications in photovoltaics Q , high frequency transistors , coherent 
phonon pulses in the description of transient thermal transport p etc. 

Time-dependent charge transport has been investigated in a number of 
studies, using different techniques. The Green-Keldysh formalism has been 
applied to transport in mesoscopic systems having external time-dependent 
voltages or barriers j^. The Lippmann-Schwinger equation was em¬ 
ployed in the context of wave packets propagation in quantum wires with 
magnetic fields [l^. More recently the time-dependent wavepacket diffusion 
(TDWPD) method was employed as an approximation to the exact stochastic 
Schrodinger equation (SSE) method H- 

In the description of the time-dependent evolution of the wavefunctions 
an essential ingredient is represented by the transparent boundary conditions 
nrocs) 0 . Similar problematics is found in beam propagation in optics 
(l6l |. As the system under direct numerical investigation is finite, the TBCs 
are required in order to ensure the wave propagation over the boundaries. 
Otherwise, the waves may be partly or completely reflected back into the 
region of interest, bringing a limitation regarding the maximum time scales 
for which the process can be investigated. These are related to the size of 
the scattering region, which, in turn, determines the computational cost. 

In this context we employ the R-matrix method as an efficient approach 
to obtain the stationary scattering functions. The formalism was developed 
by Wigner and Eisenbud 0 in the field of nuclear physics and later was em¬ 
ployed to obtain the transport properties of mesoscopic devices 0,00- 
It has been further applied to describe coherent charge transport in nan- 
otransistor s 0 . 22 . 23], thermopower in quantum wires spin dependent 
transport [2^ and to investigate the effects of graded distribution of scat¬ 


tering centers on ballistic transport |26] and cha^ localization in dendritic 
interfaces relevant for photovoltaic applications (27|. Here we use the ad¬ 
vantages of this approach as a basis for solving the time-dependent problem 
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with TBCs. 

Furthermore, an accurate description of many electronic devices of prac¬ 
tical interest usually require a quantum approach at small length scales (e.g. 
atomistic, effective mass models) and a classical description at macroscopic 
scales (e.g. drift-diffusion type models), usually combined in hybrid transport 
models, such as the drift diffusion model (QDD) B , the quantum corrected 
drift diffusion model (QCDD) or the Schrodinger-Poisson-Drift-Diffusion 
model (SPDD) js^. In this context of transport model hierarchies we discuss 
the possibility of using the developed time-dependent framework for a local 
characterization of the photo-current and charge separation that occurs in 
the vicinity of the nanostructured interfaces of photovoltaic devices. 

The paper is organized as follows. In the next section the generic model 
system is indicated and the general R-matrix formalism is presented, point¬ 
ing out the efficiency of the method in constructing the stationary scattering 
wavefunctions for a relatively large set of total energies. The time dependent 
solutions are then determined using the eigenvectors of the open quantum 
system. The computational method is discussed in detail. In the follow¬ 
ing section the charge transport across nanostructured interfaces relevant 
for photovoltaic applications is analyzed. This includes the time-dependent 
description of the charge separation and photo-current near the interface. 

2. Model and Computational method 

2.1. Coherent scattering model 

The model system, depicted in Fig. [U consists of a central finite region, 
rio) connected by semi-infinite leads, which describes a multi-terminal 
device. This framework is typically used in the Landauer-Biittiker formalism 
and can also employed in atomistic transport calculations ji^. Carriers are 
injected from each terminal and elastically scattered in the central region flo- 
The potential energy I/(r G Oq) = Fho in central region is arbitrary, while 
the potentials along the lead direction in the Vlg domains are constant, i.e. 
R(rGa) = Vh, = I/^(r^) + V;". 

2.2. The R-matrix method 

We consider the stationary Schrodinger equation 

'HT(r) = F;T(r), (1) 
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Figure 1: Two-terminal scattering model system: scattering region flo connected by leads, 
represented by fig domains. The interfaces between flo and fig are denoted by T^. A 
reference system is indicated in each lead (zg, r^). 


with H = /{2m*) A -|-t/(r), subject to scattering boundary conditions, 

i.e. the particles are coming from one of the leads 0*, commonly termed as 
asymptotic boundary conditions, since the electrons are coming from Zg —?• 
oo. The solutions of Eq. (IT]) are the scattering functions \I/(r) corresponding 
to the total energy E, which may vary continuously. 

The wavefunctions inside the leads can be generically written as: 

^g{feQs-,E) = exp {-ik^Zg)^^{rj) 

i 

+ ( 2 ) 


where ki, = are the wavevectors along the transport direction 

in each channel u. The composite index p = (s,i) denotes the channel i from 
lead s. The complex coefficients and are the amplitudes of the 
incoming and outgoing wavefunctions for each channel, with respect to the 
scattering region flo- The energies E'/_ correspond to the transverse modes $ 1 . 
inside the leads, which are solutions of the transversal Schrodinger equation 
(r^ G fig)- 




2m* 


A_l;s +V (r G fig 


= El<S>,(ri 


( 3 ) 


To obtain the solution of the scattering problem, one solves first the 
auxiliary Wigner-Eisenbud (W-E) problem, which is defined on the scattering 
region fio as the same Schrodinger equation as in Eq. ([T|) , with new boundary 
conditions on the interfaces Es! 


V-Xiir G Oo) = Eixi{r G Oq), (4) 
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( 5 ) 


dxi 

dz.. 


= 0 . 


-ir. 


One may further write the scattering functions in Oq as linear combinations 
of the W-E eigenfunctions xi as: 


eflo;E) = '^ai{E)xi{r eQo). ( 6 ) 

I 

By imposing the continuity relations on each lead - scattering region interface 
for the scattering functions and the current: 


'['(eir. 
1 cl^'(r) 


m* dzg 


1 


m* dzg 


( 7 ) 

( 8 ) 


one obtains the relation between the incoming and outgoing coefficients 
and In a compact form it may be written as where S is 


the scattering matrix. 

The S-matrix can be conveniently expressed in terms of an R-matrix 



1 

• 

1_ 

-1 

a 

s = - 

1-Rk 


1 + —Rk 


m* 


m* 


where 


with 


^ixix'{E) 


^ {xi)uix’i)u' 

2 E — El 
1=0 ^ 


ixi)i^ = f drs^u{rj)xi{r G U). 
Jr^ 


( 10 ) 

( 11 ) 


The k-matrix is diagonal, 

The scattering functions can be determined for each energy using the 
scattering S-matrix. Assuming the particle is incident on channel u, i.e. from 
lead s and having the transversal mode i, one may write the wavefunction 
inside the lead s" as: 


T^(r G fls";E) = —= W [exp {-ik^iZs') 5^^^ 

vStt , 

v' 

exp {ik^,Zs')]^u'{ri>) 5s's 


( 12 ) 
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while in the scattering region we have: 


e no]E) 




'y ^ (1 ~ G Ei) 


(13) 


with 


and 


G ^0) E) 

7?.(r G r2o) £ 


= / d^,7^(rG^o,r'G^,;£;)$,(r,^) 


r2o; E) 


^ Xiir ^ ^ ^o) 


(14) 

(15) 


Depending on the total energy E we distinguish between open (prop¬ 
agating) channels with real ky for < E^ and closed (non-propagating) 
channels with imaginary ky for E^ > E. 

The transmission functions can be determined for each pair of propagating 
modes and each energy E from the unitary matrix S = = 

The total lead-to-lead transmission can be calculated as: = 

where the summation is performed only over the open channels. 


2.3. Time-dependent problem 

Using the R-matrix method the scattering wavefnnctions can be calcu¬ 
lated for an energy E^ which may vary continuously. Being solutions of the 
stationary Schrodinger equation ([T]) the functions E) form an orthogo¬ 

nal system: 

{^,{r;E)\<if,,{r;E')) = 5,,>5{E-E'). (16) 

In practical calculations, we discretize the energy axis and use a finite set of 
energies {Ek}k- 

We assume the initial quantum state of a particle is given by d'(r,t = 
0) = d^o(r)) typically) although not restricted to, a wavepacket inside the 
scattering region. For a time independent scattering potential, the time 
evolution is described by: 


^{r,t) = Fife) exp 



(17) 


where Ck^u = (Tj^(r; £';i^)|To(r)) are the expansion coefficients of the initial 
wavefnnction in the orthonormated basis made of the scattering functions. 
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(c) 


Figure 2: (a) Zig-zag and (b) mesoporous type interfaces, inside the scattering region fio- 
The black (left hand side) regions correspond to the reference potential, while the yellow 
(right hand side) regions correspond to the band offsets, (c) Band offsets for electrons 
(Vo) and holes (Fq')- 

2.Jj.. Computational method 

The main advantage of the R-matrix method relies on the fact that total 
computational cost is divided in two parts: in the first part, the solution of 
the Wigner-Eisenbud problem {(x/, E/)}/, the transversal eigenvalue problem 
in each lead and the calculation of the overlap integrals {xi)v 

do not depend on the total energy E; in the second part, the R-matrix, the 
S-matrix and optionally transmission functions and scattering wavefunctions 
are determined for each total energy E. Usually, for moderate number of 
channels and a large enongh basis set, the Hamiltonian diagonalization per¬ 
formed in order to obtain the eigensystem {{xi,Ei)}i represents the largest 
computational cost. This step is performed only once, while the energy 
dependent calculations require a much smaller computational effort. The 
advantage of this method becomes particnlarly important when a relatively 
large nnmber of transmissions and wavefunctions need to be determined. In 
the following we detail the computational steps in our approach. 

The rank of the Hamiltonian matrix in a three-dimensional model {d = 
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3) is Nb = Nbx x Nby X Nbz, where Nbx, Nby, Nbz are the numbers of ba¬ 
sis functions corresponding to each spatial dimension. The basis functions 
4’ijkiy) = X (f>yj{y) X (f>zk{z) dehned in fig should obey the vanish¬ 

ing normal derivative imposed by the condition ([HI) at the planar interfaces 
with the leads fig. Although specihc geometries require different basis sets, a 
generic model consists of identifying Qq with a parallelipipedic box of volume 
Lx X Ly X Lz- Assuming the origin of the coordinate system is in the middle 
of flo, we may choose (pxiix) = l/\^L^cos{i7:{x + Lx/2)/Lx), with i = 1,2,... 
and similarly for (pyjiy) and (pzki^)- By a one-time diagonalization, of typi¬ 
cal cost 0{N^), the Wigner-Eisenbud functions and energies are determined. 
The transversal modes are found by solving Eq. (|31) in each lead fig, 

which are (d — l)-dimensional problems and therefore require a significantly 
lower computational cost, compared to the Wigner-Eisenbud problem. 

In the following one needs to set the maximum number of channels Ug 
considered in each lead, the total number of channels N^h = Xls being the 
rank of the R- and S-matrices, which are rigorously infinite matrices. The 
number of channels Ug should be large enough to include at least the open 
channels and an appropriate number of closed ones, to ensure the convergence 
of the calculations, in particular of the transmission functions. At this point 
one can calculate explicitly the integrals {xi)v from Eq. (fTTI) . One should 
note that all these steps are energy independent. 

In the second step, the R-matrix elements for a particular total energy E 
are found by simple summations, of total cost 0{Nb x As the k-matrix 

is diagonal, finding the matrices [l±^Rk] from Eq. ([HI) is achieved by an 
insignificant cost of 0{Nch)- The most demanding part is the computation 
of the inverse of the dense matrix [1—^Rk], which may be performed by 
standard LAPACK routines with scaling of The subsequent 

matrix multiplication in Eq. (jHI) and the similarity transform used to obtain 
the S matrix do not exceed the computational cost of matrix inversion. These 
steps are performed for a number of Ne total energies, so that the total cost 
is 0{Ne X A^ch)- -^ch ^ Nb, which is a typical situation, the computational 
burden mostly falls on solving the Wigner-Eisenbud problem. One should 
note that for an accurate representation of the transversal modes, 

the number of channels considered in the calculations should be significantly 
smaller than the number of basis elements used for diagonalizing the leads 
Hamiltonian given by Eq. (|3[). The total cost of finding the transmission 
functions is 0{N^ + Ne x N^^i) one can easily see that the cost of energy 
dependent part becomes comparable to the Wigner-Eisenbud problem only if 


the transmissions are computed for a large number of total energies, making 
the procedure effective in providing rapid varying transmission functions. 

One should point out that the R-matrix method provides the set of scat¬ 
tering functions Ek)}u,k of the open quantum system, which are con¬ 

tinuous over the entire interval, leads and scattering region. One advan¬ 
tage of our method is the possibility of including potential barriers at the 
lead/scattering region interfaces, of arbitrary height, and therefore to tune 
the contacts to the leads from fully transparent for electrons (or reflection¬ 
less) with no barrier, to almost opaque with high barriers. An even bigger 
advantage is the fully transparent case itself, which practically describes an 
infinite system with a finite Hamiltonian matrix in the scattering region and 
asymptotic wave functions outside. Alternatively, a large system can in prin¬ 
ciple be described as a finite, closed system, in a basis of wave functions 
vanishing at the boundaries. But for practical calculations the size of the 
basis must increase with the size of the system, eventually becoming pro¬ 
hibitive. Therefore, for given physical parameters and simulation time our 
method is more efficient. While the solution of the Wigner-Eisenbud prob¬ 
lem represents usually the most demanding part, the presence of the leads 
generates only a minimal additional computational effort. 

3. Application to nanostructured interfaces 

3.1. Physical systems 

Nanostructured inhomogeneous interfaces increase the light scattering 
and the device active area in photovoltaic applications, enhancing the gen¬ 
eration of photoexcited carriers. The size of the inhomogeneities should be 
comparable with the light wavelength, i.e. hundreds of nanometers. Differ¬ 
ent techniques have been used to achieve this goal, e.g. by using nanowire 
arrays or mesoporous materials. The latter have been increasingly used for 
boosting the solar cell efficiencies. As concrete examples, one can mention 
thin films of nearly spherical Ti02 or ZnO aggregates with typical sizes up to 
a micron. Furthermore, the aggregates possess inner structure, being com¬ 
posed of nanocrystallites, with typical sizes of ^15nm (^. Given the quite 
different length scales involved in the physical structures of interest and the 
operation conditions, hybrid (quantum-classical) transport models should be 
used: at a few nanometer length scales, e.g. at the nanocrystallite scale, 
a quantum mechanical transport model should be best suited, while at the 
aggregate level and larger scales a drift-diffusion model should be employed. 
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Figure 3: Zig-zag interface: scattering wavefunctions (absolute value square) for (a) lowest 
transverval mode for total energies E = 0.015 eV (upper), 0.15 eV (lower) and (b) the 8th 
mode for E = 0.09,0.15 eV. The bandoffset potential is Vq = 0.1 eV. The images include 
the scattering region and portions of the leads as described in the text. The scattering 
region corresponds to the interval (—50,50) nm on the horizontal direction. 


We consider here two different types of nanostructured interfaces, namely 
the zig-zag and the mesoporous interface, as indicated in Fig. El The first 
interface presents sharp regular variations, while the second one mimics a 
prototypical mesoporous layer. Flat-band condition is assumed for each of 
the two materials in the bulk, while at the interface there is a band offset 
Vo- Both types of structures introduce a local confinement potential for 
electrons and holes. From the holes perspective, the potentials in Fig. El are 
simply interchanged, provided the bandgap is the same in the two materials. 
Generally, for different bandgaps in the two materials, Egi and Eg 2 , the band 
offset for holes is Vq = Egi — Eg 2 + Vq. Therefore, in a non-interacting one- 
particle picture, we may describe within the same formalism the independent 
propagation of electrons and holes across the nanostructured interface 

3.2. Scattering wavefunctions 

We first analyze the solution of the stationary Schrodinger equation ([I]), 
subject to scattering boundary conditions for the two types of interfaces. 

Typical scattering wavefunctions for electrons are depicted in Figs. Eland 
m for the zig-zag and the mesoporous interface, for two different transver¬ 
sal modes, 1st and 8th, and different total energies and band offsets. The 
wavefunctions are shown for the combined system of leads (fig) and scatter¬ 
ing region (Gq). We consider a rectangular scattering region defined by the 
lengths Lx = Ly = 100 nm as shown in Fig. El and additional portions of 
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Figure 4: Mesoporous interface: scattering wavefunctions (absolute value square) for (a) 
1st and Sth modes with = 0.18 eV (upper) and Vq = 0.5 eV (lower); (b) the same modes 
and total energy, but for a lower bandoffset, Vb = 0.1 eV. 
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Figure 5; Scattering functions (absolute value square) corresponding to holes incoming 
from right, for the zig-zag (a) and mesoporous (b) interfaces. The Sth mode is represented 
at E = 0.3 eV, for a band offset Vq = 0.5 eV. 




the leads of length Lg = 50 nm are also depicted. At energies lower than 
the band offsets, there is total reflection, since the Vq potential step extends 
in the right lead to infinity. At higher energies, transmitted wavefunctions 
are observed in the right lead, together with partial reflections into the left 
lead. As one can see from Fig.IHlthe regular shaped zig-zag interface produces 
uniform patterns in the scattering functions and therefore in the charge local¬ 
ization. On the other hand, in the case of a disordered mesoporous interface, 
as one may expect, asymmetries due to randomness are present. The shapes 
of the scattering potential can be correlated to the maxima of localization, 
which indicates the feasibility of our approach in the context of describing 
fine-grain nanostructured interfaces. We also refer the reader to Ref. j27[ 
where sharp dendritic interfaces were analyzed. 

In the effective mass model, the hole perspective is mirrored with respect 
to the electron picture. Figure El shows the scattering wavefunctions for holes, 
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(a) (b) (c) (d) 

Figure 6; Snapshots of the wavepacket scattering for the zig-zag interface (a,b) and for the 
mesoporous interface (c,d), starting with to = 0 in steps of 0.04 ps from top to bottom. 
For (a,c) the bandoffset is Vo = 0.3 eV, leading to a complete reflection of the wavepacket, 
while for (b,d) Vb = 0.1 eV and partly transmitted electron waves are visible. For a better 
visibility, the colormap has been rescaled in each plot. 


which are incoming from the right hand side. For the chosen parameters, 
E = 0.3 eV and Vq = 0.5 eV, the holes are completely reflected and the two 
localization probability maps become complementary to the ones obtained 
for electrons: for the zig-zag interface, see Fig. IS^b) (upper part) and for the 
mesoporous interface, see Fig. 111(a) (lower part). 

3.3. Wavepacket propagation 

We assume wavepackets are generated near the nanostructured interface, 
corresponding to either electrons or holes. The initial wavepacket is described 
by a plane wave modulated by a Gaussian function 

•exp(-ikor), (18) 

where the function 0(r; Tq, Rc) is 1 for |r — ro| < Rc and 0 for |r — ro| > Rc, 
imposing strictly the vanishing condition for \l/o(r) beyond a cutting radius 
Rc and vl is a normalization constant. By using the 0 function one can 


To(r) = A ■ 0(r; Tq, Rc) ■ exp 


12 
















avoid introducing extra energy in the wavepacket as it may partly overlap 
with high potential energy regions. 

We consider the wavepacket is incident from the left, initially located at 
To = ( 3 :^ 0 , Vo), where xq = —25 nm, r/o = —25 nm for the zig-zag interface and 
Xo = 5 nm, = —25 nm for the mesoporous interface, i.e. it is placed in the 
middle of the largest circular nanoparticle. The wavevector ko = (lOyr/La,, 0) 
was chosen along the x direction, while a = Lx/IQ and the effective mass 
is m* = 0.0655mo (GaAs). The cutting radius is Rc = 10 nm. With these 
parameters the kinetic energy of the particle is Eq = 0.14 eV. A uniform grid 
with of 400 total energies Ej, in the range [0,0.3] eV was considered. 

Snapshots of wavepacket evolution are presented in Fig. [HI for the zig¬ 
zag interface (a,b) and the mesoporous interface (c,d). Two band offsets are 
considered, W = 0.1 eV and Vq = 0.3 eV, rendering a partially transparent 
interface and a practically opaque one, respectively. The shapes of the in¬ 
terfaces introduce different consequences in the wavepacket scattering. For 
example, the zig-zag interface tends to scatter significantly the wavepacket 
in the transversal direction, decreasing the kinetic energy on the x direction 
such that the propagation along the leads is slow. By contrast, in the case 
of the mesoporous interface, the largest nanoparticle is the primary scatterer 
and the scattered wave is more uniform. In this case the partially transmit¬ 
ted and reflected waves are moving away from the interface at a higher speed 
compared to the zig-zag interface, but some probability remains confined for 
a longer time at the mesoporous interface. 

Although the system considered here is two-dimensional, the method can 
be extended to the more computationally demanding three-dimensional case. 
The qualitative picture regarding charge localization and transport is however 
qualitatively very similar, as it was pointed out in Refs, j^, 22 . 


3.4- Charge propagation near the nano structured interfaee 

Based on our time-dependent approach with transparent boundary con¬ 
ditions we can calculate the evolution of the charge distribution and the 
photo-current. Wavepackets, which correspond to the electron/hole pairs 
generated at the interface, are scattered elastically during a time interval 
which is a typical average coherence time Tc- Our aim here is not to fully 
describe the current collected at the contacts as it may depend on other 
parameters, like e.g. recombination rates, but to describe the electron-hole 
separation at the nanostructured interface. Beyond Tc the carriers may suffer 
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Figure 7: Schematics of the photo-current model: the interface region Aq and the outer 
regions Ai and A 2 . Wavepackets are generated in Aq and the charge separation is 
evaluated for each region at a typical coherence time Tc- 

inelastic scattering processes and, at even larger time scales, they diffuse to¬ 
wards the contacts moving through the bulk of the two materials, a process 
which may be further described by a classical drift-diffusion model. 

In order to quantify the charge separation, we isolate the interface region 
Aq between positions Xi and X 2 . The outer regions Ai and A 2 , correspond 
to the two materials which make up the junction, as depicted in Fig. [71 
These three domains may possibly correspond to the Qq) regions which 
previously defined the scattering problem. We calculate the probability of 
finding the particle in the A^ domains, using the continuity relation for 
the probability current 


Qs{rc)= [ dt f f dr |T(r, t = rjp, (19) 

where j(r,t) = l/mRe[T*(r, t)(—zhV)T(r, t)] is the probability current. The 
probabilities Qs correspond to the amounts of charge separated by the inter¬ 
face, which are found in A* domains at Tc- Furthermore, one may define the 
current flowing into each A* region as Ig = dQg/dt. 

The time evolution of the charge separation is plotted in Fig. (Hlfor the 
zig-zag and mesoporous interfaces for the cases shown in Fig. El^b) and[n](d), 
in comparison with an ideally flat interface positioned at a: = 0. We define 
the regions A^ by choosing Xi and X 2 according to Fig. [7] : Xi = —40 nm, 
X 2 = 10 nm for the flat and zig-zag interfaces and Xi = —10 nm, X 2 = 40 nm 
for the mesoporous interface. The interface regions have a length of 50 nm 
along the x direction for all three structures and include the wave packet, 
which has the initial coordinates (—25, —25) nm in the case of the flat and 
zig-zag interface and (5, —25) nm in the case of the mesoporous interface. 
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Figure 8; Integrated probability flux Qg for the zig-zag (a) and mesoporous (b) interfaces. 
Qa is calculated for the regions AO (black), A1 (red), A2 (green). The dashed lines 
correspond to a flat interface. The insets show the ratio between Qo of zigzag/mesoporous 
interface and the same quantity corresponding to the flat interface, as an indicator of 
trapped charge in the Aq region. The photo-currents for each region A^ are indicated in 
(c) and (d). 

The wave packet is therefore in both cases initially placed mid-way between 
El and E 2 . 

In the following we discuss comparatively the particularities introduced 
by each type of interface. The chosen band offset of Vq = 0.1 eV allows partly 
transmitted waves. For both nanostructured interfaces we obtain, in some 
respect, a similar behavior for the Qs time dependency. In the interface region 
Qo decreases, while Qi and Q 2 increase in time and tend to saturate. The 
increase of Qi is delayed compared with Q 2 as the wave packet moves from left 
to right and gets scattered at the interface. By comparison, keeping the same 
parameters, the zig-zag interface introduces larger delays compared with the 
flat one. This interface transfers part of the kinetic energy, initially on the x 
direction, onto the y direction, leading to a slower advance speed. In contrast, 
the scattering across the mesoporous interface overall resembles better the 
data obtained for the flat interface, with one difference: Qq decays slightly 
slower. The wave is trapped at the interface for longer times, as one can also 
see from the insets of Fig. |H1 and also from sequence in Figs. IH^c) and (d), 
compared with the zig-zag interface. In addition, the increase of Qi and Q 2 is 
smoother compared to both flat and zig-zag interfaces. The lower plots show 
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the photo-currents in each region A*: the transmitted and reflected waves 
are visible as current pulses of different widths. The back scattered waves 
arrive at the Ai region at later times compared to the transmitted waves 
and the width of the current pulse is accordingly more dispersed. Although 
the time-dependent quantities Qs{t) and Is{t) depend on the specific choice 
of the As domains, the observed behavior remains qualitatively the same. 

For a complete characterization of a certain interface one has to perform 
statistical averages, taking into account the variation of the wavepacket pa¬ 
rameters (ko. To, a, shape etc). Moreover, averages on different interfaces 
from the same class would be generally necessary. To investigate the charge 
localization on either side of the interface, a distribution of Tc may also be 
considered and subsequently a drift-diffusion type model may be employed 
to calculate the collected photo-current in a concrete device. Here we pre¬ 
sented the methodology for evaluating the charge separation leaving a more 
complex analysis for a future study. 

4. Conclusions 

We introduced a general framework for describing time-dependent coher¬ 
ent transport using the R-matrix formalism. Expanding the time-dependent 
wavefunctions in the basis of scattering functions obtained by solving the 
stationary problem for the open quantum system, the transparent boundary 
conditions are introduced in a natural way. The R-matrix method provides an 
efficient procedure of calculation the wavefunctions for a relatively large en¬ 
ergy set, which is essential for an accurate description of the time-dependent 
transport. The detailed steps of the computational scheme are provided. 
As applications, we consider the scattering of wavepackets across nanostruc- 
tured interfaces, which become increasingly relevant for new generation of 
nanostructured photovoltaic devices. We discuss the transient behaviors of 
electrons crossing different types of interfaces. Finally, we introduce a model 
for calculation of the photo-current and charge localization, which provides 
a tool for optimizing nanostructured interfaces. 
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